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Abstract 

When a vortex in a two-dimensional Josephson junction array is driven by a 
constant external current it may move as a particle in a viscous medium. Here 
we study the nature of this viscous motion. We model the junctions in a square 
array as resistively and capacitively shunted Josephson junctions and carry 
out numerical calculations of the current-voltage characteristics. We find that 
the current-voltage characteristics in the damped regime are well described 
by a model with a nonlinear viscous force of the form Fp = n(y)y = j^^y, 
where y is the vortex velocity, n{y) is the velocity dependent viscosity and A 
and B are constants for a fixed value of the Stewart-McCumber parameter. 
This result is found to apply also for triangular lattices in the overdamped 
regime. Further qualitative understanding of the nature of the nonlinear 
friction on the vortex motion is obtained from a graphic analysis of the mi- 
croscopic vortex dynamics in the array. The consequences of having this type 
of nonlinear friction law are discussed and compared to previous theoretical 
and experimental studies. 

PACS numbers: 74.50.+r, 74.60.Ge, 74.60.Jg, 74.70.Mq 
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I. INTRODUCTION 



Understanding the vortex motion produced by externally applied currents has been an 
important topic in the study of the transport properties of Type II superconductors in the 
Abrikosov phase. When vortices are able to move across the system they produce Faraday 
voltages that are responsible for the I-V characteristics measured experimentally. It has 
been useful to describe the transport properties of the dilute vortex phase in terms of a phe- 
nomenological single vortex equation of motion. An example of this approach is given by 
the Bardeen-Stephens equation which has been successfully applied to conventional super- 
conductors []T[. A similar approach has been attempted in the description of the transport 
properties in two-dimensional Josephson junction arrays (JJA) P-l7[. These arrays are 2D 



lattices of superconducting islands (sites) connected by Josephson junctions (bonds). The 
unit cells (plaquettes) of these lattices can be, for example, square or triangular. In the 
JJA case the vortices are represented by eddy current patterns about a plaquette. Although 
these JJA vortices differ in several important ways from their continuum counterparts, the 
question that has been addressed by several authors is: To what extent can one use a single 
macroscopic equation of motion to describe the dynamical properties of vortices in JJA? 
Further interest in this problem has come from recent experiments in underdamped arrays 
JTT|-|r3| . These arrays were found to show hysteretic features in their I-V characteristics 



that suggest that vortices behave as particles with a mass ||11|| . Furthermore, experimental 
evidence for ballistic vortex motion was reported in triangular arrays Jl3 . 

In this paper we concentrate on the friction experienced by a JJA vortex. We investigate 
this friction in detail by numerical simulation of the dynamics of an array containing one 
single vortex. The commonly adopted vortex equation of motion assumes a frictional force 
proportional to the vortex velocity. Our results show that, instead, the friction is a nonlinear 
function of the vortex velocity that decreases as the velocity increases. We propose a new 
phenomenological friction law that accounts for the numerical results. 

Here we consider the classical regime defined by Ej >> E c = e 2 /2C, where Ej is the 
Josephson coupling energy and E c the charging energy of two islands, e the electron charge, 
and C the capacitance of the junction. In this regime quantum fluctuations are neglected, 
leaving the phases 8(r) of the Ginzburg-Landau order parameter on the islands as the only 
dynamical variables. The experiments mentioned above were reported to be in this regime. 
In this case the JJA are well-modelled by the RCS J model, defined by the total bond current 
i(r,r'), between nearest neighbor sites r and r' 

i(r,r')=p c 9(r,r') 

+ 6(r, r') + sin[0(r, r') - 2nA(r, r% (1) 

plus Kirchoff's current conservation conditions at each site. Here the dots represent time 
derivatives. The three contributions to i(r, r') are the displacement, the dissipative and the 
superconducting currents, respectively. The phase difference across a junction is 6(r, r') = 
9{r) — 9{r'). The currents are expressed in units of the junction critical current J c ; time is 
measured in units of the characteristic time l/u c = Ti/ (2eR n I c ), and [3 C = {l) c /u) p ) 2 is the 
Stewart-McCumber parameter [JH3], with the plasma frequency u p defined as <Ji = 2eI c /hC, 



and R n is the junction's normal state resistance. The bond frustration variable A(r, r') is 
defined as the line integral of the vector potential A: 
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A(r,r') = j-f r 'A-dl, (2) 

00 Jr 



with the elementary quantum of flux <fio = hc/2e. 

It has been suggested that the transport properties emerging from Eq. (|l|), in the case 
that the phase configurations in the array contain one vortex, can be described by a classical 
macroscopic model in terms of a single continuous vortex coordinate y that satisfies the 
equation of motion 0-0,0 

My + i]y = i b + i d sin (2ny), (3) 

where M = n/3 c and rj = n for a square array, and M = 2tt(3 c and rj = 2n for a triangular 
array. This equation assumes that the JJA vortex can be described as a point particle with 
mass M that, driven by a (Lorentz) force if,, moves through a sinusoidal pinning potential 
and experiences a viscous damping force with constant viscosity coefficient rj. The vortex 
mass M can be calculated |3j by equating the electromagnetic energy stored in the array 
to a vortex kinetic energy ^My 2 . The value of the depinning current, id, depends on the 
underlying lattice geometry. An estimate for id in a square lattice gives id « 0.1 while for a 
triangular lattice it is only i d ~ 0.02 |T9|| . 

If one substitutes 2ny by 9 in Eq. (|3|), one obtains the equation of motion for the phase 
difference across a single Josephson junction with current bias i&, critical current i d , shunt 
resistance 2R n (square lattice) or R n (triangular lattice) and shunt capacitance C/2 (square 
lattice) or C (triangular lattice). Eq. (|3]) has been studied extensively, mostly numerically, 
and the results show different types of nontrivial behavior depending on the values of the 
parameters in the equation [|T!|. The solutions to Eq. (|3]) exhibit a critical value % = id 



above which the junction is in a nonzero voltage state. In the array case i d corresponds to 
the depinning current above which the vortex moves and a finite voltage is measured. This 
voltage is proportional to the vortex velocity and arises because the phase differences across 
the array change in time. When M — 0, the overdamped case, and for currents % < id 
the vortex is pinned to the lattice and the voltage is zero. For % > id the vortex can move 
under the action of the current and a nonzero voltage state is produced. For M ^ the 
I-V characteristics resulting from Eq. (|3]) show hysteretic behavior when the current % is 
ramped up and down past the depinning current id. If M is sufficiently large Eq. (|3|) predicts 
ballistic vortex motion, in the sense that a high-velocity vortex would continue its motion 
over many lattice constants when the driving current is switched off. Van der Zant et al. 
|T3[ reported experimental evidence of ballistic vortex motion in a region without driving 
currents in a H-shaped triangular array with (3 C = 46. The quoted (3 C value was computed 
from the normal state resistance of the junctions. At low temperatures and voltages, the 
effective /3 C , determined by the quasi-particle resistance, can be orders of magnitude larger 
13| . In contrast to this experimental result, in numerical simulations within the RCSJ 



model, no evidence for ballistic vortex motion has been found [|r4}-|TT|j ; a high velocity vortex 
in an underdamped array does not move more than one plaquette as soon as the driving 
current is switched off. Furthermore, the calculated I-V characteristics for square arrays 
show almost no hysteresis near the vortex depinning current for f3 c = 10 [|l(J, whereas Eq. 
(0) would yield a substantial hysteretic behavior. 

In trying to understand this discrepancy between experimental and theoretical studies 
an additional dissipative mechanism, arising from the coupling of the vortex to spin waves 
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or plasma oscillations, has been shown to give rise to a nonzero vortex viscosity in the 
completely under damped limit (/3 C = oo) 1 15fl . This would invalidate the model (§) in this 



limit, leading to very small mean free paths over which vortices come to rest if the driving 
current is switched off, even in highly underdamped arrays. The enhanced viscosity has also 
been measured experimentally in Ref. Nevertheless is was suggested there that the 

vortices might still move ballistically in a current-free region at low velocities. 

In this paper we carry out a systematic comparison between the results for the I-V char- 
acteristics obtained from simulations of JJA described by Eq. ([!]) and the I-V characteristics 
obtained from an equation of the form given in Eq. ([|). The analysis is carried out for a 
range of (3 C values. We find that an equation of the form of Eq. (^) is not representative of 
the JJA results. Instead we find strong quantitative evidence that an equation that yields 
a rather good fit to the JJA results is 

M(p c )y+ y = k + i d sH^y)- (4) 

Here the constants A, B and M are found to be weakly dependent functions of (5 C . This is 
the main result of this paper. We note that the linear friction law given in Eq. (|3]) has to be 
modified in a nonlinear way to account for the JJA results. This nonlinear dependence on 
the vortex velocity applies in particular to the range j3 c = up to (3 C ~ 100. This change to 
a nonlinear dissipation law raises some important questions, for example, how to introduce 
temperature effects at the phenomenological level f20fl . We will discuss other important 



consequences emerging from this nonlinear viscosity law later in the paper. 

The outline of the paper is as follows. In Section II we discuss the calculational algorithm 
used to compute the I-V characteristics from Eq. (|1]). In Section III we present the bulk 
of our results for the I-V characteristics together with the fitting analyses that lead to the 
result given in Eq. (^). In this section we also discuss the microscopic aspects of the vortex 
motion in the array by analyzing the current distributions of the vortex as a function of time. 
Section IV contains our conclusions together with a comparison to previous experimental 
and theoretical work. 



II. CALCULATIONAL APPROACH 

In this paper we are interested in calculating the dynamical response of an array of 
Josephson junctions driven by a constant d.c. current. The set of nonlinear dynamical 
equations of motion given in Eq. ([!]) can be efficiently integrated using a fast Fourier 



transform algorithm |21 



In our simulations we use a square lattice (with a lattice constant set equal to unity) 
with periodic boundary conditions (pbc) along the y-direction while the current is fed in and 
taken out along the x-direction (see Fig. 1). Hence a vortex tends to move in the ^-direction. 
The total number of plaquettes along the x- and y-directions are denoted by N x and N y , 
respectively, whereas the total number of sites is L x x L y , with L x = N x + 1 and L y = N y . 
Most of the results presented in this paper correspond to systems with L x = L y . However, 
these results do not change significantly when considering systems with L x < L y . In fact, 
if L x is not smaller than 8, the weak finite size effects encountered in our calculations are 
mainly governed by the vortex motion along the y-direction. 
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The vorticity n(R) of a plaquette R can be defined as (see for instance the Appendix of 
Ref. |U) 

27rn(Jl) = 2nf + J2 (^( r > r ') ~ ^A(r, r')). (5) 

V{R) 

Here V(R) denotes an anticlockwise sum around the plaquette -R and the gauge invariant 
phase difference 6(r, r')—27rA(r, r') is taken between — 7r and +tt. The frustration parameter 
/ measures the average flux piercing a plaquette, measured in units of 4>q. Physically, vortices 
in JJA can be seen as eddy currents in the current flow pattern. If there is only one vortex 
in the array, then there is one plaquette, say Rq, with vorticity n(R ) = 1 while all other 
plaquettes have zero vorticity. We will call Rq the core of the vortex. In a JJA with pbc 
the phase configurations corresponding to a single vortex in the middle column of the array 
cannot be written down as easily as in an array with free boundaries. We construct a 
single vortex configuration with a method used previously in Ref. ||. It allows for a direct 
calculation of the phase configuration in terms of the vorticities n(R) G {—1,0, 1} once a 
gauge choice for the A(r, r') and a choice for one of the phases 9{r) has been made. Here 
we are mainly interested in understanding the one- vortex dynamics and thus we concentrate 
on this case throughout the paper. 

We take the frustration equal to / = l/N x N y so that the single vortex we introduce in the 
middle column of the array has a current pattern symmetric around that column. The single 
vortex equation of motion proposed in Eq. (|3J) or Eq. (^) describes a continuous motion in 
the ^-direction, whereas the location of a vortex as determined from the phase configurations 
is discrete and undetermined within the vortex core. Therefore, when making comparisons 
between the vortex velocity, defined in terms of the microscopic phases, to that obtained 
from the coarse grained vortex variable y, we need to compare time averaged quantities. 
The vortex velocity is directly related to the time average of the voltage V(t) across the 
array in the x-direction, where V(t) is defined as 

V(t) = L j2^[e(L x -l,y)-9(0,y)], (6) 

y=0 al 

according to the Josephson relation. Time is again measured in units of 1/uj c , and V(t) is 
measured in units of R n I c . Each time the vortex has travelled over a distance N y the total 
phase difference across the array has changed by 2n and therefore the vortex velocity v is 
given by 

V = W 
In 

where V =< V(t) >. 

Another quantity of interest in describing the vortex dynamics is the current vorticity 
around a plaquette defined as 

C(R,t)= £i(r,r',t). (8) 

V(R) 

An important difference between the vorticity n(R,t) and C(R,t) is that the former is an 
integer, while the latter is a continuous function describing the vortex as an eddy current 



5 



pattern extending outside the vortex core. One interesting quantity to look at is the "center 
of mass" of the current vorticities 

~ _ E R R V C(R) _ E Ry R y D(Ry) 
Af Af 

with 

D(R y )=Y,C(R). (10) 

Rx 

The normalization factor Af has an unusual form depending on the pbc assumed in our 
calculations. It is determined by the requirement that Y v has to change one unit if the current 
vorticity configuration is shifted one plaquette. For two current vorticity configurations 
C(R) and C'(R) shifted by exactly one plaquette with respect to each other, we get 

Af = Y,RyD\R y )-Y,RyD{R y ) 

Ry Ry 

= Y J D{Ry)-L y D{L y -l). (11) 

By taking large enough lattice sizes, one can ensure that the quantity Af is essentially 
constant for a range of positions of the vortex core in the middle of the coordinate system. 
In this region Y v show steps with integer height magnitudes. 

III. RESULTS 

In this section we present the evidence we have found that allows us to conclude that a 
vortex in a JJA moves with a nonlinear viscosity law, at least in the overdamped {j3 c = 0) 
to damped (/3 C < 100) regime. We start by considering the overdamped case, in which there 
are no shunt capacitors, and therefore there is no spin-wave dissipation channel. Next we 
will discuss the results for values of /3 C up to 100. 

A. Nonlinear viscosity in the /3 C = case 

Fig. 2 shows a typical I-V characteristic computed for a 32x32 array with one vortex 
(/ = l/N x N y = 1/992) and with (3 C = 0. In order to compare with the result from Eq. ([3|), 
we have plotted the vortex velocity as defined in Eq. (|7|) versus the current. The results 
for the I-V characteristics for larger lattices are practically the same, as we will discuss in 
more detail below. For currents if, < 0.10, the vortex is pinned by the lattice and thus the 
measured voltage is zero. Above if, pa 0.10, the vortex is depinned from the lattice by the 
current and its motion gives rise to a nonzero voltage. Up to approximately if, = 0.97, the 
time-averaged vortex velocity gradually increases and so does the measured voltage across 
the array. This is the regime where the phenomenological equation of motion must apply 
and thus we call the range if, G [0, 0.97) the vortex regime. At if, > 0.97 we enter a current 
regime in which eventually all individual junctions in the current direction perform phase 
slips. 
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(9) 



In Fig. 2 we also show, as a continuous line, the result for the I-V characteristic assuming 
the validity of Eq. (Q), with the identification of the parameters as deduced in Refs. 
In the M = case considered here the known analytic expression for the I-V characteristic 
is given by, 



V = ^- = 2*Z*. (12) 

/of*/ * 

In Fig. 2 we observe that the vortex equation of motion seriously underestimates the time- 
averaged voltage almost everywhere in the vortex regime. There is quantitative agreement 
with the calculated results only for bias currents very close to the depinning current. More 
interestingly, we observe a qualitative difference between the two curves. Whereas the 
viscosity rj in Eq. (|3]) is constant the simulations show an effective viscosity which decreases 
with increasing vortex velocity. This leads us to propose as a model a vortex equation of 
motion of the form (j3 c = 0), 

V(y)y = ib + idsin(2ny), (13) 

with a velocity-dependent viscosity rj = rj(y). We have found that the functional form for 
rj(y) that fits our results in the vortex regime quite well reads, 

V (y) = A/(l + By). (14) 

where the sign of y is taken positive and A and B are parameters determined by fitting the 
array results to this form. As in the constant viscosity case we can analytically evaluate the 
result for the I-V characteristics yielding, 



il 



A-BJi 



(15) 



The top curve in Fig. 3 shows a fit using Eq. fllS]) in the current range ij, = 0.10 to % = 0.80, 
to the P = I-V characteristic obtained from the simulations of a 32 x 32 array. The values 
for the parameters are A = 2.67 and B = 1.80. To indicate the error bars of these values, 
we mention that, if we fit the form ( |i~5|) to the simulation results in the range iy, = 0.10 to 
ib = 0.70, the fitted values for A and B are approximately 0.5% and 2.0% larger, respectively. 

To check on a possible size dependence of these results we carried out the same analysis 
for lattices of sizes 64 x 64 and 128 x 128. We find that the I-V characteristics leads to 
essentially the same results for currents in the range ib G (0.10, 0.80). Representative results 
from the fits for 16 x 16 and 32 x 32 arrays are given in the inset of Fig. 3. The results for 
the two larger sizes are essentially indistinguishable, within their error bars, from the ones 
shown for 32 x 32 arrays. 

We conclude at this stage that the (3 C = results are rather well described by the 
phenomenological Eqs. ( |T3"D and flT4]). In the inset of Fig. 5 we show the comparison between 



the friction force Fd in the constant viscosity model (Eq. (Q)) and the one proposed here 
(Eq. (Hp) using the /3 C = values for the parameters A and B derived from our fits. 
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B. Nonlinear viscosity in the (3 C ^ case 



We move now to discuss the changes that occur when (3 C 7^ 0. We solve the phenomeno- 
logical vortex equation of motion, Eq. numerically to compare to the results obtained 
from solving the JJA equations. In Fig. 3 we show the results for the I-V characteristics of 
a 32 x 32 array in the vortex regime for /3 C = 0, 1.2, 2, 3 and 25. We also show the fits to 
the array results using the form given in Eq. (|l^). Nonlinear behaviour persists for values 
of f3 c up to 50. The fits for small f3 c are of the same quality as for (5 C = 0. Note that we take 
M(f3 c ) = in these fits, although Eq. (|3|) suggests a non-zero mass M as soon as (3 C 7^ 0. 
Including a mass M as a parameter, as in equation (4), does not result in a better fit of the 
nonlinearity in the vortex regime. The choice M(/3 C ) = is corroborated by the fact that 
we do not find any measurable hysteresis near the depinning current in the simulated I-V 
characteristics for Stewart-McCumber parameters even up to (3 C ~ 35, in agreement with 
the results for (3 C = 10 reported previously in Ref. [To]. The puzzling conclusion is then that, 



even when the microscopic equations of motion have a "mass" term, the phenomenological 
vortex equation of motion behaves as if the vortex mass is zero or very small. In the inset 
of Fig. 3 we show the j3 c dependence of the parameters A and B for (3 C up to 4. We note 
that A increases slightly with f3 c while B decreases slightly . This trend indicates that the 
viscosity becomes "more" linear as f3 c increases. This trend can be understood as being a 
consequence of the "spin-wave friction" mechanism that sets in at (3 C > and leads to an 
enhancement of the linear viscosity. An approximate estimate of the linear viscosity in this 



regime presented in Ref. [15] led to a rise roughly proportional to y(3 c . The same trend was 



found experimentally and a semi-quantitative explanation of the results was given in the 



second reference of 11 



For lattices of size 8x8 and larger one needs (3 C values of the order of 100 to detect 
small hysteresis loops near the depinning current. For these lattice sizes, no hysteresis is 
measured up to f3 c = 35 in our , using a current grid as small as 3 x 10~ 6 . We have found, 
however, that for a small array of size 4 x 4, a very small hysteresis loop is visible in this /3 C 
regime which resembles in shape the ones obtained using Eq. (|3]). In Fig. 4 we show these 
hysteresis loops for (3 C values between 7 and 13. Note that all the I-V characteristics have 
the same depinning current while ramping the current up whereas they have different zero 
voltage intercepts when lowering the currents. 

In Fig. 5 we show the vortex regimes of the I-V characteristics for f3 c values up to 100 
on a 32 x 32 lattice. Here we only show the result from ramping up the current, thereby 
omitting the small hysteresis loops below the depinning currents for (3 C = 50 and 100. For 
ftc > 7, we find sharp jumps in the voltage. The vortex regime ends at these jumps, which are 
believed to be due to switching of rows of junctions to the resistive state. This row-switching 



behavior has been seen before in experiments |illJT2| , |22[| and in simulations fi4 Hl7| , [23f . In 



this figure we note a crossover from a nonlinear to a linear viscosity regime as j3 c increases. 
At (3 C = 100 the vortex regime of the I-V characteristic is nearly linear for if, > 0.25 and can 
be extrapolated through the origin (the step-like structure of the I-V characteristics in the 
upper half of the vortex regime corresponds to interference of the vortex with its periodic 
image, as was explained in Ref. ]l5j , and will disappear if we consider a system with larger 
L y ). For (3 C = 50 a similar extrapolation does not intersect the origin, so the friction is still 
nonlinear. We note that the range of applicability of the nonlinear viscosity model given in 



8 



Eq. (^]) covers some of the (3 C values reported in the experiments in Refs. 
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C. Nonlinear viscosity in triangular arrays 

All of the calculations described above were performed in square lattices. Recently Yu 



and Stroud carried out calculations of the I-V characteristics in triangular arrays [pL7| . We 
have reanalyzed their results in light of our nonlinear viscosity model given in Eq. (||). In 
Fig. 6 we show the corresponding fit to their (3 C = results for the I-V characteristic of an 
8x8 lattice with the current biased in the [101] direction. The fits to the parameters A and 
B yield the results 7.67 and 2.47, respectively. From these numbers we conclude that the 
magnitude of the viscosity is roughly 2.8 times as large as in the square lattice case, while 
the nonlinearity parameter B/A is smaller. We note that a theoretical prediction of a factor 
of two for A between the square and triangular lattices was made in Ref. (9). 



D. Microscopic vortex motion 

It is interesting to directly study the time evolution of the vortex motion across the 
array in order to further understand the nature of the nonlinear viscosity suggested by the 
phenomenological macroscopic model given in Eq. We concentrate here on the f3 c = 
case. 

In Fig. 7(a) we show current vorticity distributions of a vortex for % = 0.11, slightly 
above the depinning current, at different times. We observe that the vortex motion has 
essentially two time scales, a slow and a fast one. In the slow regime the vortex does not 
move much while it gets deformed by the applied current. Subsequently the vortex moves 
fast until it gets stuck again and the stretching process repeats itself. This type of stick-slip- 
like motion is reflected in the nonlinearity of the viscosity. The decrease of the viscosity with 
increasing velocity is analogous to the behaviour of the kinetic friction coefficient between 
dry surfaces in the stick-slip phase, which is likely to be generic for frictional dynamics at 
low speeds |[24| . As shown in Fig. 7(b) the qualitative motion of the vortex remains the 



same for % = 0.6, although the quantitative values for the slow and fast times have become 
smaller. 

In Fig. 8 we show results for the time-dependent voltage V(t) (full line) and the normal- 
ized center of mass vortex velocity (dashed line) defined as 

5(f) = jY v , (16) 

where Y v is given in Eq. (Bp, for three values of %. In order to ensure a constant value for the 
normalization factor M in a wide range of vortex positions we choose a lattice of size 8 x 64. 
We observe that in both the V(t)- and the w(t)-curves, the amplitude of the oscillations 
around the average value decreases with increasing the bias currents. If we interpret the 
quantity v(t) as a coarse-grained vortex velocity, the physical meaning of this result is that 
the pinning force decreases when the vortex velocity increases. To check this interpretation 
we extract the pinning barrier from the simulations, by measuring the variation in the array 
energy given by 
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E = cos[0(r,r') - 2nA(r,r')}. 

<r,r'> 



(17) 



We find that the pinning barrier shows a similar decrease as a function of the (time-averaged) 
vortex velocity. Under the assumption that the pinning force in the array is proportional 
to the pinning energy barrier, we conclude that indeed the amplitude of the pinning force 
decreases when the vortex velocity increases, in accordance with the interpretation of the 
quantity v. We can qualitatively relate this result to the current vorticity snapshots shown 
in Fig. 7: for larger currents the vortex moves faster, the current vorticity spreads out over 
more plaquettes and the pinning at the core plaquette becomes less effective. 



IV. CONCLUSIONS AND COMPARISON TO PREVIOUS WORK 

In this paper we have proposed a phenomenological vortex equation of motion that fits 
well the I-V characteristics obtained from solving the full set of JJA microscopic dynamical 
equations. The main difference with previous studies is that our proposed equation of motion 
has a nonlinear velocity dependent viscosity that decreases as the velocity increases. The 
validity of this description covers the range from overdamped to damped regimes as defined 
by the Stewart-McCumber parameter, and it applies to square as well as to triangular 
lattices. We also have found that for (3 C < 35 the the I-V characteristics indicate that the 
vortex moves as if its inertial mass is zero, or at most very small (no hysteresis at depinning). 
As j3 c increases the nonlinearity of the viscosity slowly decreases at the same time that the 
linear term slowly increases. We will now discuss the above results in the light of previous 
experimental and theoretical studies. 

Experimentally, evidence for the nonlinear viscosity can be seen in the I-V characteristics 



reported in Refs. |TT]jr^| . However, the I-V characteristics results measured in Ref . |TT| for an 
almost overdamped triangular array do not show evidence for a nonlinear viscosity, whereas 
in our simulations it is in the overdamped case that the nonlinearity is dominant (see Fig. 
5). 

On the theoretical side, Eckern and Sonin || derived a general vortex equation of motion 
in the continuum limit. In the adiabatic, or small vortex velocity limit, this equation reduces 
to the model of Eq. (§) without the sinusoidal pinning force. This equation of motion is 
believed to take into account the spin- wave friction occurring when f3 c ^ 0, as found in Ref. 



15| . Here we will focus on the f3 c = case, in which the equation of motion also provides 
corrections to Eq. (|3p beyond the adiabatic limit. Taking a constant vortex velocity v, for 
a constant current bias, for (5 C = the vortex equation of motion reduces to 



f 7,2 p -k/V^ 

"/ a '#tt^t 2 ™" (18) 

Here the integral in fc-space is over the two-dimensional plane. The vortex velocity v is 
taken along the ^/-direction. The exponential in the integrand provides a smooth cutoff for 
large k. An alternative cutoff used in Ref. consists of replacing the exponential in the 
integrand by the two-dimensional Heaviside function 

e(\k x \-7r)e(\k y \-ir). (19) 
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Note that in Eq. ([181) we cannot add a sinusoidal pinning force, as this would be inconsis- 
tent with the constant-velocity assumption. However, the inclusion of the pinning poten- 
tial, being most important in producing a finite depinning current, would barely affect the 
higher-velocity part of the I-V characteristics, where the nonlinearity in Eq. fll8|) is most 
pronounced. 

In Fig. 9 we show the I-V characteristics computed from Eq. fllTf) , with the exponential 
fc-cutoff, and for the cutoff given in QI5p. Both curves intersect the origin, because the lattice 
pinning potential is absent in fll8l) . We note that the inclusion of non-adiabatic effects in 
this equation of motion gives rise to a viscosity that decreases with increasing velocity. 
Although there is an improvement in the higher-velocity part of the I-V characteristic when 
compared to the linear viscosity vortex equation of motion, Eq. (||), the predictions of the 
continuum model still deviate qualitatively from the full (lattice) calculations. We note that 
the higher-velocity component of the I-V characteristic depends crucially on the choice of 
high-momentum cutoff in Eq. flTB|). 

This work has been motivated in part by the issue of ballistic vortex motion. The phe- 
nomenological vortex equation of motion presented in this paper attributes a mass M((3 C ) = 
to the vortex in square arrays in the regime of j3 c < 35. This is a consequence of the absence 
of hysteresis in the I-V characteristic in this regime, also reported in Ref. |TBJ. This means 



that the electromagnetic energy stored in the shunt capacitors does not represent a kinetic 
energy for the vortex in this regime, at least not in the way suggested by the model Eq. 
(0). This detracts from the idea behind the possibility of ballistic vortex motion in JJA 
described by the RCSJ model. 

In a separate argument, the enhancement of the viscosity with increasing j3 c , leads to 
very small path lengths over which a vortex with high initial velocity looses its assumed 
kinetic energy. The enhancement of the viscosity was also measured experimentally [JTT 



It has been explained in Refs. [|Pli|l5| in terms of an additional friction mechanism due to 



coupling of the vortex to plasma oscillations. It was also suggested in Ref. || that this 
coupling would not prevent ballistic vortex motion in a small velocity window in triangular 
arrays. Recent simulations |HJ of triangular arrays did not show such a velocity window 
in the parameter range considered (0 < f3 c < 1000). However, one may need much larger 
values of [3 C to possibly see ballistic vortex motion Hl3~| . 

For the discrepancy between the results of the experiment of Van der Zant et al. [13j 
and that of the simulations based on the RCSJ model, one possible explanation suggested 
recently in Ref. [ID] involves the discreteness of the charges in the array. The clarification of 
this problem needs further experimental and theoretical study. With regard to the nonlinear 
vortex viscosity found in our work, establishing a direct connection between the microscopic 
and the phenomenological description represents a difficult problem for future study. 
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FIGURES 



FIG. 1. Array geometry used in the simulations, illustrated with an 8 x 8 array. Junctions are 
denoted as crossed bonds. In the y-direction periodic boundaries are imposed, while the current 
bias is applied along the x-direction. 

FIG. 2. (3 C = I-V characteristics, plotted as average vortex velocity versus normalized bias 
current. The dashed line gives the results from simulations of a 32 x 32 array with one vortex. 
The simulation shows a vortex viscosity that decreases with increasing vortex velocity. The full 
line was obtained from the model vortex equation of motion (3). 

FIG. 3. Simulation results for the I-V characteristics as in Fig. 2 for different values of (5 C 
(circles). The full lines are the fits to the I-V characteristics using Eq. (15) for a 32 x 32 lattice. For 
clarity of presentation the origin of successive (3 C values is offset to the right by 0.1 unit. The inset 
shows the values of the fitted parameters A and B as a function of (3 C and array size. Diamonds 
(A) and triangles (B) correspond to 32 x 32 array whereas circles (A) and squares (B) to 16 x 16. 

FIG. 4. Hysteresis loops in the simulated I-V characteristics for a 4 x 4 array with one vortex, 
for different /3 C values. Note the smallness of the current scale. First the current is swept up to 
% = 0.1265, slightly above the depinning current, and subsequently it is swept down until the 
vortex is retrapped by the lattice. The depinning current is % = 0.126 for all f3 c values shown. 
For j3 c = 7 the current at which the vortex is retrapped is (on this scale) equal to the depinning 
current. For higher values of (3 C this current is increasingly lower. For 8x8 and larger lattices the 
hysteresis for these (3 C values disappears. 

FIG. 5. I-V characteristics (vortex velocity versus bias current) from simulations of a 32 x 32 
array with one vortex, for different values of (5 C . From top (marked with T') to bottom ('2'), 
P c = 0, 1.2, 3, 7, 15, 50, 100, respectively. Note that the higher /3 C , the smaller the vortex regime. In 
the inset the nonlinear friction force Fr,{y) = Ay /{I + By) (full line) is shown as a function of the 
vortex velocity y, for the f3 c = parameter values A = 2.67 and B = 1.80. The dashed line is the 
friction force Fr,{y) = vry in the constant-viscosity model (3). 

FIG. 6. Simulation results from Ref. [17] for an 8 x 8 triangular array with current along the 
[101] direction, [3 C = (circles). The continuous line is the fit using Eq. (15). The fit parameters 
obtained are A = 7.67 and B = 2.47. 
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FIG. 7. Snapshots of a smooth interpolation of one vortex current vorticity distribution in an 
8x8 sublattice of a 16 x 16 array, for two different current values: (a) % = 0.11, (b) % = 0.60. 
Different gray scales represent different levels of current vorticity. In the first snapshots (labelled 
as '0'), the black dot is the middle point of the vortex. The dashed line (fixed in time) is a guide to 
the eye. In (a) the time interval between two snapshots is At = 10 (in units of l/u c ). The vortex 
moves over one plaquette in approximately t = 7 At. In (b) the time interval between successive 
snapshots is At = 0.625 in units of 1/lo c , and here the period of the motion is approximately 
t = 4.5A*. 

FIG. 8. Rescaled voltage (full lines) and vortex center of mass velocity (dashed lines) versus 
time, in an 8x64 array, for three different bias currents. The amplitude of the oscillatory component 
in both quantities decreases with increasing %. 

FIG. 9. Comparison of the f3 c = average vortex velocity versus normalized bias current 
obtained from a simulation of a 32 x 32 array with one vortex (dashed line) with the one from Eq. 
(18) (curve (a)). Curve (b) is obtained from Eq. (18) by replacing the smooth high-momentum 
cutoff in (18) by the sharp cutoff (19). Curve (c) is obtained from the model vortex equation of 
motion (3). 
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